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We develop a method based on the collisionless Boltzmann equation to calculate the gravitational 
clustering of relic neutrinos in realistic cosmological models dominated by cold dark matter (CDM) 
and the cosmological constant. This method can be used to estimate the phase-space distribution of 
any light particles in CDM halos. We find that neutrinos with masses > 0.3 eV cluster appreciably 
in dark matter halos above the galactic size. The resulting neutrino overdensity above the cosmic 
mean neutrino density increases with both the neutrino mass and the halo mass, ranging from ~ 10 
for 0.3 eV neutrinos in ~ 10 13 Af© halos to ~ 1500 for 1.8 eV neutrinos in ~ 10 15 M Q halos. We 
examine the implications of neutrino clustering for the Z-burst model of ultra high energy cosmic 
rays (UHECR), which interprets the observed events at E > 4 x 10 19 eV as decay products of 
Z-bosons from the resonant scattering between relic and high energy neutrinos and anti-neutrinos. 
We estimate the UHECR energy spectrum for various neutrino masses towards five of the most 
massive clusters in the local universe (within 100 Mpc): Virgo, Perseus-Pisces, Hydra, Centaurus, 
and Coma. The UHECR flux in the Z-burst model is expected to be significantly higher towards 
these clusters if m„ > 0.3 eV and nearly isotropic otherwise. 



I. INTRODUCTION 

The nature of cosmic rays above the Greisen-Zatsepin- 
Kuzmin (GZK) cutoff [1] at - 4 x 10 19 eV is an un- 
solved problem in ultra high energy cosmic ray (UHECR) 
physics [2]. These events have been reported by the 
Akemo Giant Air Shower Array (AGASA) [3], Fly's Eye 
[4], Havera Park [5], HiRes [6], and Yakutsk [7] collab- 
orations. Interactions with the cosmic background pho- 
tons 7 cm & via photoproduction of pions (jyy cm b ~~ * P + 
Nit, TVYcmf, — > n + Ntt), photopair production (pj cm b —> 
pe + e~ , jjcmb ~> e+e )i and inverse Compton scatter- 
ing (e ± 7 cm6 ->• e ± 7, pj cm b -* Pi) at high energies [8] 
constrain a ~ 10 20 eV cosmic ray to a few Mpc for the 
characteristic lengths of either charged cosmic rays or 
neutrons and photons. More specifically, the attenuation 
length of protons above the GZK cutoff is ~ 50 Mpc. 
The lack of known processes to accelerate cosmic rays in 
small Galactic objects makes the Galactic origin of these 
ultra high energy particles unfeasible [9] . Novel powerful 
acceleration mechanisms for light nuclei are required if 
these energetic particles are produced in nearby galax- 
ies [10]. Exotic particles and dynamics have also been 
suggested, but these come with their own difficulties [2]. 

One proposed explanation for the UHECRs is the Z- 
burst model, which tries to solve the puzzle without in- 
voking new physics beyond the standard model of particle 
physics except for neutrino masses. Several recent exper- 
iments [11-14] have found evidence for non-zero neutrino 
mass. The Z-burst model hinges on the fact that ul- 
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tra high energy neutrinos (and anti-neutrinos) produced 
at cosmological distances can reach the GZK zone unat- 
tenuated. Their resonant annihilation on the relic anti- 
neutrinos (and neutrinos) produces Z bosons, about 70% 
of which decay into hadrons within ~ 10 -25 sec. The final 
state has fifteen pions and 1.35 baryon-antibaryon pairs 
on average [15], where the fifteen pions decay into thirty 
high energy photons. The Z boson is highly boosted 
(~ 10 10 ) [16], resulting in a highly collimated beam with 
a half angle of ~ 10~ 10 . This and the fact that the effect 
of magnetic fields at such high energies is negligible [17] 
ensure a high probability for the protons and photons to 
reach the observer if the Z-burst occurs in the direction 
of the Earth. The Z-burst model has been discussed in 
detail in many papers [18-24]. The resulting cosmic ray 
flux has been shown to depend strongly on the density of 
the relic neutrinos [15, 16, 20, 25, 26], but the neutrino 
density in these calculations has been taken to be either 
the constant relic density from the big bang or some ad 
hoc value. 

In this paper we perform a detailed calculation of the 
neutrino clustering in the local universe using realistic 
cosmological models and apply the results to the Z-burst 
model for UHECRs. Since the current constraints from 
cosmological observations and laboratory experiments in- 
dicate that the neutrino masses are small (<^ 2 eV; see 
Sec II) and the CDM dominates the dark matter den- 
sity (f2 c dm 3> ^i/), we do not expect the clustering of 
neutrinos to affect significantly that of the CDM. As a 
result, it is not essential to use full scale, time consuming 
iV-body simulations. Instead, we solve the collisionless 
Boltzmann equation for the neutrino phase space distri- 
bution in a background potential given by the universal 
profile of CDM halos reported in recent high resolution 
simulations [27]. The Boltzmann equation is then linear 
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in the neutrino density contrast and has tractable inte- 
gral solutions. The advantage of this method over the 
conventional iV-body simulations is that we can obtain 
the neutrino density profile much below the resolution 
scale (~ 50 kpc) of large cosmological simulations by us- 
ing as an input the CDM potential determined from much 
higher resolution simulations of individual halos. More- 
over, the computation time required for our approach is 
negligible compared with numerical simulations, thereby 
allowing us to explore a large parameter space of neutrino 
masses and dark matter halo masses. 

In Sec II the relevant Boltzmann equation and the in- 
tegral solutions are derived. In Sec III results for the 
clustering of neutrinos for different neutrino masses and 
CDM halos are presented and compared with physical 
arguments based on neutrino free streaming and the 
Tremaine-Gunn constraint [28]. The resulting neutrino 
density profiles are also compared with earlier TV-body 
simulations [29], which show good agreement. In Sec IV 
the neutrino overdensity calculation is incorporated in 
the Z-burst model for UHECRs, where we make realistic 
predictions for the UHECR energy spectrum for different 
neutrino masses. We estimate the level of anisotropy in 
the UHECR flux by examining lines of sight towards five 
of the most massive clusters (Virgo, Perseus-Pisces, Hy- 
dra, Centaurus, and Coma) in the local universe (within 
100 Mpc) where neutrinos are likely to be most clustered. 



II. BOLTZMANN EQUATION FOR NEUTRINO 
CLUSTERING IN CDM HALOS 

In this section we develop an approach based on the 
collisionless Boltzmann equation to study how massive 
neutrinos cluster gravitationally in realistic cosmological 
models. We start by noting that the median velocity of 
unclustered background neutrinos of mass m„ (in eV) at 
redshift z is 
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This implies that light neutrinos (m„ < 2 eV) do not 
accrete significantly onto CDM protoclusters until z ~ 3 
because the neutrino thermal velocities are greater than 
the velocity dispersion of a typical cluster or superclus- 
tcr. We are then faced with the more tractable problem 
of how neutrinos cluster in the potential well of an exist- 
ing CDM halo. Our approach is to use the collisionless 
Boltzmann equation for the neutrino phase space density 
/ and follow its evolution in a background CDM poten- 
tial given by the approximate universal profile obtained 
in high resolution simulations of individual halos [27]. 
Note that the CDM potential is time-dependent in gen- 
eral. Earlier work has used the Boltzmann approach to 
study how neutrinos cluster around point masses in the 
context of cosmic string seeded galaxy formation [30, 31]. 
This method allows us to calculate the neutrino density 
profiles in the inner part of the cluster (<J 10 kpc) that 
can not be resolved by large cosmological simulations. 



This will be seen to be important in the Z-burst model 
where a significant contribution to the cosmic ray flux 
comes from the inner regions of the halo. 

In the Newtonian approximation and in physical co- 
ordinates, the collisionless Boltzmann equation takes the 
form 
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dt 



r- V r f + p- V p / = 0. 
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Rewriting it in conformal time dr = dt/a and in comov- 
ing position X = r/a and momentum q = ap m v dr, 
we obtain 
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where the Newtonian gravitational potential $ obeys 
p = — m„V<I>. At the time of decoupling the neu- 
trino phase space density is given by the thermal Fermi- 
Dirac distribution fo(q) °c {e q l Tv > a + 1) _1 where T v $ = 
(4/11)1/3^^ = i 6 7 6 x W -i eV is the temperature of 
the neutrino background today. Gravitational cluster- 
ing distorts the spatially uniform /o, so we write the full 
distribution as 

f(x,q,r) = f (q) + f 1 (x,q,r). (4) 

The gravitational potential can also be written as 

&(x,t) =$ (x,t) + $ 1 (x,t), (5) 

where $o is related to the mean background comoving 
density p = p cdm + p u by V x $ = ^-Gpa 2 X, and $i 
is determined by the density contrast of both CDM and 
neutrinos in the halo: 



V^$i = 4irGa 2 (6p cdm + Sp„) 
Eq. (3) then becomes 



(6) 



ldf 



— t\ 1 ■ V x /i - ra.V^i • V q fo 

a or 



-m v V x $i-V q fi=0, (7) 

where we have used the Friedmann equation d = 
— ^-Gpa, which gives aaX + V z <l>o = 0. We note that 
eq. (7) is the full Boltzmann equation and no approxima- 
tion has been made thus far. It is generally a nonlinear 
equation in /i where 4>i is related to f\. For our problem, 
however, <I>i is mostly determined by the CDM whose po- 
tential has a well known, pre-determined form. Eq. (7) 
is then linear in f\ and much easier to solve. 

Furthermore, eq. (7) has a simple integral solution if 
the fourth term is neglected. For example, in earlier cal- 
culations that examined neutrino clustering onto point 
masses seeded by cosmic strings, this term was dropped 
in order to simplify the calculation [30, 31]. We will 
also drop this term, but we justify this approach in two 
ways. First, we note that dropping this term requires 
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Vg/i < V 9 /o and not fi < /q. The former is generally 
a less restrictive condition and can be satisfied even if fi 
is much larger than f because on dimensional grounds, 
we have 



Vg/i ^ fi/cr v ^ Sp„/a* 
V 9 /o /o/w pV/« 4 



5, 



(8) 



where a v is the velocity dispersion of the gravitational 
potential $i and v is the median neutrino thermal veloc- 
ity in cq. (1). Since only neutrinos with v a v are cold 
enough to fall into the gravitational wells, we expect the 
ratio V g /i/V 9 /o to be much smaller than /i//o, thereby 
making it easier to justify the dropping of the fourth 
term in eq. (7). (For example, V 9 /i/V,j/o ^ 0.2 for <~ 1 
eV neutrinos in ~ 10 14 M Q halos.) Eq. (8) further indi- 
cates that the neutrino overdensity is larger than f\ / /o by 
(a v /v) 3 , a large factor in highly clustered regions. This 
explains qualitatively the large overdensity found in our 
numerical results to be presented in Sec III. In the next 
section we also provide further justification for ignoring 
the Vg/i term by comparing our results with the neu- 
trino density profiles of two halos obtained from earlier 
numerical simulations. 

Following [30] , we convert eq. (7) into an ordinary dif- 
ferential equation by going into Fourier space and using 
a new time variable dr) = dr/a: 

where /i and p are the Fourier transforms of /i and p. 
The solution is 
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where we have taken the initial neutrino phase space to 
be homogeneous, i.e., /i(r?o) = and /(jy ) = jo- 
The comoving neutrino number density is given by 

h v {k,ri) - n v = A j d 3 qh{k,q,rj), (11) 

which can be obtained from eq. (10) after integration by 
parts in q and integration over angles: 



n v {k,rj) - n v = 



L 



no 



^ Mk X~X ] l mv] (i2) 



where n v w 112 cm -3 is the cosmic mean comoving num- 
ber density of one species of light neutrinos and antineu- 
trinos, and h p is the Planck constant. Eq. (12) is the 
main equation that we will solve in this paper. It is 
a Volterra integral equation of the second kind that has 
the form f(t) = f* dsK(t, s)f(s)+g(t) and can be solved 



with the trapezoidal rule. It describes how neutrinos of 
a given mass m v cluster around a realistic CDM halo as 
a function of time. The density of the halo at a given 
time, p, should generally be the sum of the CDM and 
the neutrino components, but as we verify below, using 
the CDM potential alone is a very good approximation 
for the cosmological models of current interest. 



III. RESULTS FOR NEUTRINO CLUSTERING 
IN CDM HALOS 

In this section we present results for n v computed from 
cq. (12). We choose to integrate eq. (12) from z = 3 to 0; 
the results differ by only about 5% if the initial redshift 
is pushed back to 5 because the neutrinos do not cluster 
appreciably at such early times as discussed above. We 
also need to specify the cosmological models and neu- 
trino masses. For the cosmological parameters, we use 
the currently favored critically-flat model with matter 
density fl m — fi c dm + Cl v = 0.35, cosmological constant 
f^A = 0.65, and Hubble parameter h = 0.7. Variations in 
these parameters at 10 to 20% level are not expected to 
alter our neutrino results significantly since the effect on 
the halo potential <&i is small. For the neutrino masses, 
we consider four different models in which three mod- 
els assume three degenerate massive species with masses 
0.6, 0.3, and 0.15 eV respectively, and one model with 
one massive species with mass 1.8 eV. The correspond- 
ing density parameter in neutrinos is f2„ ss 0.04, 0.02, 0.01 
and 0.04, respectively, all much smaller than f2 c d m - This 
range of neutrino masses is chosen to span the current 
cosmological and laboratory limits. The most recent cos- 
mological constraint comes from the galaxy clustering 
power spectrum of 2-Degree-Field Redshift Survey, which 
places an upper limit of 1.8 eV on the sum of the neutrino 
masses [32]. The Supcr-Kamiokande experiment [11] pro- 
vides strong evidence for oscillations between neutrino 
species with a mass difference of Sm 2 = (1 — 8) x 10~ 3 
eV, giving a minimum mass of ~ 0.07 eV if the neutrino 
masses are hierarchical. The choice of three degenerate 
neutrino masses is based on indications that if any of the 
mass eigenvalues is above 0.1 eV then all three masses 
are above 0.1 eV and almost degenerate [33]. 

Although the final neutrino density profile will depend 
strongly on the CDM profile, we do not expect the re- 
verse to hold: CDM density will not be much affected by 
the clustering of neutrinos because all the models con- 
sidered in this paper have small Cl v /Cl c dm, so the CDM 
dominates the gravitational potential of a dark matter 
halo. We are therefore justified in using the universal 
CDM profile determined from the high resolution pure 
CDM simulations [27] as the input: 



PcdmM = 



p5rl 



r(r + r s ) 2 ' 



(13) 



where 8 and r s are given in terms of the concentration 
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parameter [34, 35] 



c = 



M 



1 + z \ 1.5 x lO^-iM© 
200c 3 



-0.13 



3 [ln(l + c) - c/(l + c)] ' 



(14) 



1.63 x 10~ 5 



M 



ft" 1 Mr 



1/3 



/i _1 Mpc. 



We evaluate the mass density p on the right hand side of 
eq. (12) exactly by adding the CDM density given above 
and the neutrino density computed from previous time 
steps. We did find that approximating p with the CDM 
profile alone (i.e. ignoring the neutrino contribution to 
the total potential) changes our results by no more than 
10% for the light neutrino masses and cosmological pa- 
rameters considered in this paper. We also tested the 
simplifying assumption made in [30] , which allowed them 
to reduce the integrals in eq. (12) analytically to a sin- 
gle integral by using [e q l Tv - a + 1) _1 = Ae~ q l Tv > a , i.e., by 
assuming a Maxwell-Boltzmann instead of Fermi-Dirac 
distribution. We found this simplification to cause a ~ 
20 % difference at small scales, so we do not use this 
approximation here. 

Before presenting our results for the realistic models 
above, we first conduct a comparison study by checking 
our results for the neutrino density profile against those 
in the numerical simulations of Ref. [29], which inves- 
tigated the clustering of CDM and neutrinos in two flat 
cosmological models: f2 c( jm = 0.8 and Q v = 0.2 with two 
species of 2.3 eV neutrinos, and f2 c dm = 0.7 and fi„ = 0.3 
with one species of 7 eV neutrinos. Both models assumed 
h = 0.5. These models are no longer consistent with cur- 
rent observations, but the simulation results provide a 
useful tool for us to test the validity of our Boltzmann 
approach. For a fair comparison, we use the CDM halo 
profile found in [29] as an input: 



PcdmM = 



c 



r(r + R) a ' 
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where R and a are 0.3 Mpc and 1.5 for the f2„ = 0.2 
model and 0.11 Mpc and 1.1 for the Vt v = 0.3 model. We 
note that their outer profile is shallower than eq. (13) 
from the pure CDM simulations [27]. In Ref. [29] this 
difference was attributed to the differing spectral index 
n for the matter power spectrum of the neutrino models 
at the mass scale of the simulated halos (w 1.3 x 10 15 
M Q ): n ~ -1.36 for fi„=0.2 and n ~ -1.53 for 0^=0.3. 
This argument, however, appears inconsistent with the 
near universal nature of the halo density profile reported 
in Ref. [27]. A better understanding for the origin of halo 
profiles should help resolve this issue. 

Fig. 1 compares the ratio of the neutrino and CDM 
density profiles from our approach vs. the two simulated 
halos in Ref. [29] . For both halos we have used the same 
cosmological parameters, CDM profiles, and halo mass 
in our Boltzmann approach as in their simulations. We 



find a good agreement between the two methods for the 
outer parts of the cluster; whereas our results are lower 
by about 50% in the inner parts. This discrepancy may 
be due to our neglecting the fourth term in eq. (7) , or due 
to statistical fluctuations in the substructure in their sim- 
ulations. (Only two simulated halos are presented in Ref. 
[29].) The Boltzmann approach used here also allows us 
to explore how neutrinos respond to different CDM po- 
tentials. As an illustration of this, we show in Fig. 1 how 
the neutrino profile changes when the input CDM profile 
is changed from eq. (15) to the higher resolution profile 
of eq. (13). We conclude from Fig. 1 that we can obtain a 
reasonable estimate for neutrino clustering using eq. (12) 
instead of full scale iV-body simulations. 

We now turn to our results for the realistic cosmolog- 
ical models and neutrino masses given at the beginning 
of this section. The four panels in Fig. 2 show the neu- 
trino overdensity n v jn v computed from eq. (12) for four 
models of neutrino masses. The more massive neutrinos 
clearly cluster more because of their lower thermal veloc- 
ities. Within each panel, the four curves illustrate how 
n v increases with halo masses from 10 12 to 10 15 M Q as 
a result of the deeper halo gravitational potentials. The 
growth of n v in the inner parts of the halo, where it is 
almost independent of r, is illustrated in Fig. 3 for 0.7 
and 0.4 eV neutrinos in 10 15 ,10 14 , and 1O 13 M halos. 
Most of the clustering is seen to occur at low redshifts. 

Unlike the CDM density that continues to rise towards 
the inner part of a halo as p oc 1/r, all curves in Fig. 2 
show that n v flattens out at some radius. Similar features 
were also seen in Ref. [31] for neutrinos clustered around 
cosmic strings. This relative suppression in the neu- 
trino vs. CDM overdensity reflects neutrino free stream- 
ing, which dampens and retards perturbation growth on 
small length scales due to phase mixing. The neutrino 
damping scale, Rd, can be characterized by the length 
scale above which neutrinos behave like the CDM. The 
standard method to solve the Boltzmann equation for 
a fluid with pressure invloves the transformation of the 
Boltzmann equation into an infinite hierarchy of velocity 
moment equations [36], where the lowest three moments 
with I = 0,1, and 2 correspond to the density, velocity, 
and shear of the fluid. The choice of the truncation of the 
hierarchy depends on the physical properties of the fluid 
and the length scales. For CDM, for example, all modes 
above I — 1 are zero. The parameter Rd gives the scale 
above which the Boltzmann hierarchy for neutrinos can 
be truncated at I = 1 (as for the CDM), and below which 
more Z-modes must be included to compute the neutrino 
damping effect accurately. This scale is given by [37] 



y/l + [a(r)/a„ r ] 2 



(16) 



where a nr ~ ST^ o/m^ is the expansion factor at which 
the neutrinos become non-relativistic. For the cosmo- 
logical models considered in this paper, the scale is 
Rd(z = 0) w 5.8/to„(cV) Mpc. From Fig. 1, we in- 
deed find the ratio p v j ' p c dm to be about the cosmic mean 
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value Qv/flcdm at scales above Rd and to decrease grad- 
ually at smaller radii (top panels), with a final flattening 
in the neutrino overdensity at r ~ 0.1 Rd (bottom pan- 
els). Fig. 2 shows that the radius at which 8 V flattens 
out depends weakly on the mass of the CDM halo. It 
occurs at smaller radii for less massive CDM halos pri- 
marily because the lower mass halos provide shallower 
potential wells. The damping scale Rd of eq. (16) is 
to be contrasted with the neutrino free streaming dis- 
tance, which is typically defined as the comoving dis- 
tance traversed from the time of neutrino decoupling to 
a nr : X fs = J^ r dT'/a{T r ) ~ 600/m„(eV) Mpc [38]. The 
distance A f s reflects the global streaming motion of neu- 
trinos but not the local clustering properties of neutrinos 
around CDM after they become non-relativistic. 

Another way to understand the results in Figures 1 and 
2 is to compare the thermal velocities of neutrinos with 
the velocity dispersions of the CDM halos: neutrinos can 
cluster significantly only if their mean thermal velocity 
in eq. (1) is below the typical velocity of the host CDM 
halo. Fig. 4 compares these two characteristic velocities 
for a range of neutrino masses and halo masses. Since the 
NFW profile specifies only the spatial and not the veloc- 
ity distribution of the CDM halo particles, two velocity 
ellipsoids are shown for comparison: isotropic, which is 
appropriate near the center of the halo, and the more 
radial distribution (3 = 1 — Vt/ V r = 0.5, which is appro- 
priate for the outer regions. Fig. 4 illustrates that < 0.15 
eV neutrinos are too hot to be captured significantly by 
< 10 14 M Q halos, while the more massive neutrinos can 
fall into progressively lower mass halos, a result consis- 
tent with that shown in Fig. 2. 

Our results for neutrino clustering can be compared 
with the Trcmaine-Gunn bound [28], which gives an up- 
per limit on the neutrino density in the core of a halo 
based on the argument that the maximum coarse grained 
phase space density can not exceed the maximum ini- 
tial phase space density due to phase mixing. Their re- 
sults are not directly applicable to our problem because 
in their derivation, neutrinos are assumed as the sole 
constituent of dark matter, and the coarse grained neu- 
trino distribution is assumed to be Maxwell-Boltzmann 
instead of Fermi-Dirac for computational convenience. 
More recent work [39] has extended the derivation to 
models including both CDM and neutrinos and obtained 
Pv < |2<I>| 3 / 2 m 4 /127r 4 , where $ is the gravitational po- 
tential of the system. For the NFW profile, we find 
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where r s and 6 are the CDM halo parameters given by 
eq. (14). For 1.8 eV neutrinos, for example, this formula 
gives n v /n u < 3.9 x 10 4 ,3.2 x 10 5 ,2.7 x 10 6 , and 2.3 x 
10 7 for 10 12 , 10 13 , 10 14 and 1O 15 M halos, respectively, 
at the scale radius r s . One can see that this constraint 
is satisfied by at least three orders of magnitude for all 
neutrino overdensities in Figures 1 and 2. 



IV. IMPLICATIONS FOR ULTRA HIGH 
ENERGY COSMIC RAYS 

In this section we apply the neutrino clustering results 
from Sec III to the Z-burst model for UHECRs. No pre- 
vious work on the Z-burst model has included realistic 
calculations for n v . Instead, the value of n„ has been 
chosen based on certain observational constraints [24] or 
physical arguments [16, 20] and has differed greatly from 
n v ~ (1 - 10 5 )h,y. For instance, in [24], it is inferred from 
the CDM distribution in our local universe, but the large 
smoothing scale ~ 20 Mpc assumed in the calculation 
results in n v ~ n v . In contrast, our results from Sec III 
show that n v can be ^> n v in the inner ~ 1 Mpc of CDM 
halos. In Rcfs. [16, 20], n v is approximated based on 
phase-space arguments similar to that of [28] . While this 
approach gives an upper bound on the neutrino cluster- 
ing, the actual overdensity can be significantly less, as 
we have discussed in the previous paragraph. In addi- 
tion, the neutrino clustering scale of « 5 Mpc assumed 
in [20] is much larger than what we find in our calcula- 
tions. Our method gives specific predictions for n v as a 
function of halo radius, halo mass, and neutrino mass. 

To estimate the cosmic ray flux we follow the stan- 
dard assumption in the Z-burst model that the UHECRs 
above the GZK cutoff are produced by the resonant vv 
scattering, while the lower energy events are explained by 
protons originating from a uniform distribution of extra- 
galactic sources. The latter appears consistent with the 
isotropic distribution of E < 4 x 10 19 eV events detected 
in AGASA and HiRes [40]. We compute the cosmic ray 
flux (in (eV m 2 s sr) _1 ) from the Z-burst model with [24] 



F Z (E) = I dE, 
lo 



] p I 
Jo 



dr 



/ dE Vi F Vi {E Vi , r)n Di (r) + / dE Vi F- Ut (E Di , r)n Vi (r) 

JQ JO 



xovp(s) Bt(Z — > hadrons) 



dN, 



p+n 



dE r , 



dP p (r,E p ;E) 
dE 



(18) 
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Here F Vi [E Vi , r) is the flux of ultra high energy neu- 
trinos with energy E Vi at distance r and n Vi {r) is the 
physical number density of the relic neutrinos. (The re- 
peated index i is summed over different neutrino species.) 
The particle interactions are described by the cross sec- 
tion ovp(s) for the Z-boson production at the center-of- 
mass energy s = 2m v E v , and by the branching ratio 
Br(Z -> hadrons) = 69.89 ± 0.07 % for the subsequent 
cascade of the Z-boson into hadrons [41]. The factor 
dNp+ n /dEp gives the energy distribution of the produced 
protons and neutrons. The subsequent proton propaga- 
tion is specified by P p (r, E p ; E), which gives the prob- 
ability that a proton created at distance r with energy 
E p arrives at Earth with an energy greater than E. It 
measures the amount of proton energy degradation due 
to the resonant photoproduction of pions and other pro- 
cesses discussed in Sec I. Specific values of P p has been 
calculated in [42] for the range of parameters considered 
in this paper. We do not include in our UHECR flux 
estimate the contributions from the photons produced in 
the Z-decay because experimental data suggest that less 
than 50% of the cosmic rays above 4 x 10 19 eV are pho- 
tons at the 95% confidence level [43, 44]. Typical phys- 
ical mechanisms used to explain the suppressed photon 
contributions are large universal radio background and 
sufficiently strong extragalactic magnetic fields (> 10~ 9 
Gauss) [24, 45] . The study of the effects of these param- 
eters on the UHECR flux due to the Z-burst model can 
be a subject of future work. 

A key feature of the Z-burst model is that the cross 
section a u p(s) for v 9 — ► Z° is enhanced by several orders 
of magnitude near the resonant energy in the rest frame 
of the relic neutrinos [18] 



2to„ ( 



4.2 x 10 21 eV 



lcV 



in,, 



(19) 



where Mz is the mass of the Z boson. The flux in eq. (18) 
to a good approximation therefore depends only on the 
neutrino resonant energy and not on the slope of the 
incident high energy neutrino spectrum. Eq. (18) can 
then be written as [24] 



F Z (E) = a„ p F Vi (E™ s ) \ dE, 



dr 



dP p (r,E p ;E) 



(I jq 



dE 



(20) 



where n Vi (r) is the physical number density of neutrinos 
and antineutrinos at the Z-burst site at radial distance r, 
£T, y p = 40.4 nb is the cross section for v D — » Z° averaged 
over the width of the resonance, and F Vi {E^ a ) is the in- 
cident flux of ultra high energy neutrinos at the resonant 
energy. The function Q p is the boosted momentum dis- 
tribution from hadronic Z decays and can been calculated 
from experimental data [24] . It has a fairly broad peak at 



y PS 10~ 2 and falls off approximately as y~ 7 for y > 0.5. 
The neutrino flux F Ui (El es ) remains a free parameter in 
the Z-burst calculation since no successful astrophysical 
model yet exists to explain the production of > 10 21 eV 
neutrinos [46-48] . We do not attempt to model the effect 
of source evolution in our calculations since it is again an 
unknown quantity and is easy to incorporate once its na- 
ture is known. 

We first present results for the cosmic ray spectrum 
F(E) in the Z-burst model ignoring the spatial clustering 
in the neutrinos, i.e., we assume n v = n v in eq. (20). This 
assumption underestimates the flux in the Z-burst model, 
but we include the results here for comparison since this 
is a common assumption made in several Z-burst calcula- 
tions [16, 24, 48]. Fig. 5 shows the sensitivity of E 3 F(E) 
on the neutrino masses. The flux is higher at high ener- 
gies for smaller m v because the momentum distribution 
Q p of the decay particles peaks at a higher energy for 
smaller m v . For a given m„, E 3 F(E) decreases rapidly 
at E ;> 10 21 eV because Q p falls off as ~ y~ 7 for y > 0.5. 
The integration is carried out to a maximum distance of 
Rmax — 2000 Mpc, but our results are insensitive to this 
choice as long as R m ax is sufficiently beyond the GZK 
zone of ~ 50 Mpc. 

For ease of comparison, the curves in Fig. 5 are all nor- 
malized to the same incident neutrino flux of F Vi (i?£ es ) = 
1.7 x 10~ 35 (eV m 2 s sr) _1 for each of the three neutrino 
flavors. (For the one flavor 0.07 eV model, the assumed 
flux is 3 times higher.) We do not attempt to determine 
this value by performing statistical fits to data because 
the UHECR spectrum from AGASA (square symbols) 
and HiRes (triangle symbols) disagree in both amplitude 
and shape. We do note that for models that have three 
degenerate neutrino masses of m Vi ^ 1 eV, this value 
for the neutrino flux is consistent with the existing up- 
per bound from the Goldstone Lunar Ultra-high energy 
neutrino Experiment (GLUE) [49]. The 0.07 eV model 
shown in Fig. 5, however, would need to be lowered by a 
factor of ~ 4 in order to satisfy the GLUE upper limit. A 
better understanding of systematic effects in the GLUE 
experiment is needed before their results can be used to 
rule out models. 

For comparison, the dotted curve in Fig. 5 shows the 
cosmic ray flux for protons originating from a uniform 
distribution of extragalactic sources with a constant co- 
moving density. It is computed from 



Feg(E) 



dE v 



xF p {Ep) 



R„ 



dr 



Pmax 

dPp{r,E p ;E) 



[l + z(r)} 3 



dE 



(21) 



where the unknown proton injection energy spectrum 
Fp(Ep) is typically assumed to be a power law: F p {E p ) = 
e~ 1 AE~ 13 . The shape of Feg depends on the injec- 
tion spectrum F p , but for definiteness, we have assumed 
j3 = 2.4 and A = 5.98 x 10 31 (and an upper energy cutoff 
of E p = 10 23 ), which are found to be the best fit values 
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[24] to the existing cosmic ray data that have a total ex- 
perimental exposure of e w 8 x 10 16 m 2 s sr. The GZK 
cutoff is clearly seen at <~ 4 x 10 19 eV in the dotted curve. 
The flux rises beyond ~ 4 x 10 20 eV because the photo- 
production of pions is a resonant process where the cross 
section peaks at E p ~ 2.3 x 10 20 eV [1] and decreases at 
higher energies, allowing a larger fraction of protons to 
reach us. 

The predictions for the UHECR spectrum change sig- 
nificantly when we incorporate the neutrino overdensity 
computed in Sec III. To make realistic estimates for our 
local universe, we consider five lines of sight towards 
five of the most massive nearby clusters: Virgo, Cen- 
taurus, Hydra, Perseus-Pisces, and Coma, where the 
highest overdensity of neutrinos are expected. The dis- 
tance, mass, rough angular extent, and equatorial co- 
ordinates of each of the clusters are listed in Table 1 
[50, 51]. The cluster masses are taken from http://cfa- 
www.harvard.edu/~huchra/clusters, where they are esti- 
mated from galaxy velocities and the virial theorem. Wc 
caution that these values have large error bars. The mass 
of the nearest cluster Virgo [52-55], for example, has 
been estimated to be 1.5 — 6 x 10 14 M Q based on X- 
ray emission measurements by ROSAT [52], to 1.5 x 10 15 
M0 based on the relativistic Tolman-Bondi method [53] . 
(The Tolman-Bondi model is based on analytic solutions 
to the Einstein field equations for spherically symmet- 
ric pressure free overdensities in a homogeneous universe 
[53, 56].) 

Fig. 6 shows our predictions for the cosmic ray flux (in 
(eV m 2 s sr) _1 ) towards these five lines of sight for four 
different neutrino masses. Along each line of sight, high 
energy neutrinos from extragalactic sources are assumed 
to traverse a uniform sea of background neutrinos plus 
an overdensity of background neutrinos centered at the 
location of the given cluster, where n v is computed from 
cq. (12) for the mass of the cluster. We also include in 
our calculation a local n v for the Local Group of mass 
4 x 10 12 M Q [57]. Despite the smaller mass, the proximity 
of the Local Group to us leads to non-negligible contri- 
butions to the UHECR flux: about a factor ~ 2 for the 
3 x 0.6 cV model, and up to a factor of ~ 10 at E > 10 20 
eV for the 1.8 cV model. The difference is primarily due 
to the more efficient clustering of 1.8 eV neutrinos com- 
pared to 0.6 eV neutrinos in the Local Group. 

Our main conclusion from Fig. 6 is that the flux of 
UHECRs in the Z-burst model should show significant 
anisotropy if m„ 0.3 eV, with the largest flux com- 
ing from the Virgo cluster. For m u < 0.1 eV, on the 
other hand, neutrinos are too hot to cluster appreciably 
even around the largest clusters in the universe, and the 
UHECR flux in the Z-burst model is nearly isotropic. We 
choose to plot in Fig. 6 the flux per steradian because the 
angular extent of the clusters cannot be precisely defined, 
but one can easily use the approximate angular extents 
of the clusters listed in Table 1 to estimate the expected 
anisotropy in the signal. 



V. CONCLUSION AND DISCUSSION 

We have introduced and tested a method based on the 
collisionless Boltzmann equation to calculate the gravita- 
tional clustering of massive neutrinos in CDM halos for 
realistic cosmological models. This method is valid for 
currently favored models with f2 c d m 3> &v in which the 
clustering of neutrinos is mostly determined by the exist- 
ing CDM halos while the clustering of the CDM is little 
affected by the neutrinos. One can then obtain the neu- 
trino phase space distribution by solving the collisionless 
Boltzmann equation in a background potential given by 
the universal profile of CDM halos from high resolution 
simulations. The resulting Boltzmann equation is linear 
in the neutrino density contrast and has tractable inter- 
gral solutions that require negligible computational time 
in comparison with N-body simulations. This method 
has enabled us to obtain specific predictions for the neu- 
trino overdensity as a function of halo radius, halo mass, 
and neutrino mass for a wide range of parameters. 

Our calculation shows that neutrinos with masses j> 
0.3 eV can cluster appreciably in CDM potential wells 
with masses > 10 13 M Q . The predicted neutrino over- 
density increases with both the neutrino mass and the 
halo mass, ranging from ~ 10 for 0.3 eV neutrinos in 
- 10 13 M© halos to - 1500 for 1.8 eV neutrinos in 
~ 1O 15 M . Specific predictions are plotted in Figs. 2 and 
3. Neutrino clustering has a strong impact on the Z-burst 
model that has been proposed as a possible explanation 
for the UHECR events. The predicted UHECR spec- 
trum shown in Figs. 5 and 6 depends sensitively on the 
neutrino mass and overdensity, showing distinct spectral 
features towards nearby galaxy clusters if m v > 0.3 eV. 

To illustrate the effects of neutrino mass and overden- 
sity on the UHECR spectrum, we have chosen to nor- 
malize the flux in Figs. 5 and 6 with the same value (i.e. 
F vi (Kf) = !-7 x 10~ 35 (cV m 2 s sr)" 1 for each flavor for 
the three degenerate mass models and three times higher 
for the one massive species model) instead of adjusting it 
by fitting individual spectrum to existing data. We have 
nonetheless included current data from the AGASA [3] 
and HiRes [6] experiments in Figures 5 and 6 for com- 
parison. More events are needed to discriminate the dif- 
ferent models and the directional dependence. The large 
increase in flux towards Virgo is an interesting signature 
of the Z-burst model for upcoming experiments such as 
Auger [58] and OWL [59] that will provide an angular res- 
olution of ~ 1°. Experimental limits on the anisotropy 
would in turn imply small neutrino inhomogeneities in 
the Z-burst model and can be used to place upper bounds 
on the neutrino mass. 

A useful constraint on the Z-burst model is provided 
by the Energetic Gamma Ray Experiment Telescope 
(EGRET) measurement of the GeV 7-ray background 
flux, which must not be exceeded by the high energy 
photons produced in the Z-burst models once they cas- 
cade down to the GeV energy range. The result depends 
on the assumed redshift evolution of the sources that pro- 
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TABLE I: Parameters of five nearby clusters 



Name 



Distance(Mpc) 



Mass(M ) 



Angular Radius(°) 



RA(h min) 



Dec(° 



Virgo 
Centaurus 
Hydra 
Perseus-Pisces 
Coma 



15 
43 
53 
76 
99 



7.9 x 10 14 
1.3 x 10 15 

4.6 x 10 14 
5.5 x 10 15 

1.7 x 10 15 



5 

1.5 
2 
7 

2.5 



12 29.6 
12 46.1 
10 34.5 
03 15.3 
12 57.4 



+11 49 
-41 02 
-27 16 
+41 20 
+28 15 



duce the incident high energy neutrinos, and on whether 
the sources themselves produce photons. The normaliza- 
tion of the neutrino flux cited in the previous paragraph 
rules out sources emitting a comparable flux in 7-rays 
because it leads to a conflict with the existing EGRET 
limits for the GeV 7-rays. For pure neutrino sources, cal- 
culations based on particle transport codes show that for 
neutrino masses of 0.1, 0.5, and 1 eV (ignoring neutrino 
clustering), the EGRET bound is met for a <, —3, <, 0, 
and <; 3, respectively, where the source number density 
evolves as (1 + z) a [24, 60]. When neutrino clustering is 
taken into account, our results from Fig. 2 show that the 
bound above for m v <, 0.3 eV should be unaffected since 
they do not cluster appreciably in the Local Group. For 
larger neutrinos masses, however, we expect a less strin- 
gent bound on the source evolution due to local neutrino 
clustering. To derive quantitative constraints would re- 
quire detailed transport calculations. 

The implications of the neutrino clustering results pre- 
sented in this paper extend beyond the problem of the 
UHECR spectrum. For UHECR, upcoming experimen- 
tal results may converge on a spectrum that is consis- 
tent with the GZK cutoff and would therefore not re- 
quire models such as the Z-burst. It is also likely that 
the Z-burst model is not the correct explanation for the 



UHECR events. However, the neutrino-anti-neutrino res- 
onance scattering process remains one of few ways to de- 
tect the relic neutrinos, as first suggested in Ref. [18]. 
This paper has addressed neutrino clustering, a major 
uncertainty in all studies concerning relic neutrinos. 
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FIG. 1: Neutrino clustering calculated with our Boltzmann approach (dashed) vs. numerical simulations from [29] (dotted). 
The simulation resolution is 62.5 h" 1 kpc. The upper panels show the ratio of the neutrino mass density p„ to the CDM 
density p c dm as a function of radius for two halos of 1.3 x 10 15 M Q in two cosmological models. The lower panels show p v /p v 
for the same models. The solid curves compare neutrino clustering around CDM halos with an NFW profile [27] to illustrate 
how neutrinos respond to different gravitational potentials. 
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FIG. 2: Total neutrino number density n v (r) as a function of halo radius for different neutrino masses and halo masses at 
redshift 0. The curves are all normalized to n v ~ 112 cm~ 3 for ease of comparison. The four panels (from left to right) show 
how the clustering decreases as the neutrino mass is lowered, a result of increasing neutrino thermal velocity and more effective 
free streaming. Within each panel, the curves show how n„ decreases as the halo mass is lowered from 10 15 to 10 12 M Q , a result 
of shallower gravitational wells and smaller halo velocity dispersions compared with the neutrino thermal velocity. This figure 
shows that neutrinos with m v ^ 0.15 eV cluster appreciably in M ;> 10 12 M Q halos. 
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FIG. 3: Time evolution of the neutrino overdensity in the inner parts of the halos (where n„ is independent of radius) for 
m v = 0.7 (left) and 0.4 eV (right). In each panel, three halo masses 10 15 , 10 14 , and 10 13 Mq are shown (top down). Neutrinos 
start to cluster significantly only at late times, with > 75 % of the clustering taking place between z = 1 and 0. 
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FIG. 4: Velocity dispersion of NFW halos of mass 10 15 , 10 14 10 13 , and 10 12 Mq (top down) as a function of halo radius. Two 
velocity orbits for the halo particles are shown for comparison: isotropic (solid) and (3=1 — vljv^ — 0.5 (dot-dashed). The 
horizontal lines indicate the present-day median thermal velocity for 0.15, 0.6, and 1.8 eV neutrinos. The values suggest that 
rn v ^> 0.15 eV neutrinos are cold enough to cluster gravitationally, particularly in massive halos. 
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FIG. 5: Predictions for the cosmic ray flux produced in the Z-burst model when the relic neutrino density is assumed to be the 
big bang uniform background density without gravitational clustering. Three of the four models shown assume three degenerate 
masses, each with 0.6, 0.3, and 0.15 eV (from bottom up); the fourth model assumes a single massive species with 0.07 eV. The 
cosmic ray spectrum of protons originating from a uniform extragalactic background sources is shown for comparison (dotted). 
The GZK suppression in the flux is clearly seen at E <; 4 x 10 19 eV in all spectra. The squares show the current 30 UHECR 
events from AG AS A [3]; the triangles show the HiRes events [6]. 
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FIG. 6: Predictions for the cosmic ray flux produced in the Z-burst model using realistic neutrino overdensity computed in 
this paper. The four panels compare the UHECR spectrum for the same four neutrino mass models as Fig. 5. Within each 
panel, our predictions for the spectrum towards five of the most massive clusters in the local universe are shown: Virgo (solid), 
Centaurus (dotted), Hydra (dot short-dashed), Perseus-Pisces (short dashed), and Coma (long dashed). For comparison, the 
dot- long-dashed curve shows the spectrum when neutrino clustering is ignored (i.e. the same as in Fig. 5). For ra v <; 0.3 eV, 
we predict that the Z-burst model should produce distinct spectrum towards each line of sight. For m v <^ 0.1 eV, neutrino 
clustering is insignificant and the spectrum is expected to be nearly isotropic as seen in the lower right panel. The data points 
are the same as in Fig. 5. 



